;;;;;; This note gives some basic operations of matrix and vectors


(library-directories '("/home/cl/gits/libs-for-chezscheme"))
(import (range)
        (parallelization-in-vector))

(define (vector-apply proc v)
  (apply proc (vector->list v)))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; vector
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; note: 并行计算应在最需要的地方使用,而不是处处!

;; sum elements of a vector
(define (vector-sum v)
  (vector-apply + v))
(define (vector-sum/parallel v)
  (vector-apply/parallel + v))

;; norm
(define (vector-norm v)
  (sqrt (vector-sum (vector-map * v v))))
(define (vector-norm/parallel v)
  (sqrt (vector-sum (vector-map/parallel * v v))))

;; cv
(define (vector-scalar-product c v)
  (vector-map (lambda (ele) (* c ele)) v))
(define (vector-scalar-product/parallel c v)
  (vector-map/parallel (lambda (ele) (* c ele)) v))

;;; matrix: vector of vectors, e.g. (vector (vector 1 2 3) (vector 4 5 6)) is 2X3 matrix
(define (matrix-ref A i j)
  (vector-ref (vector-ref A i) j))

;; jth column of a matrix
(define (matrix-ref-column A j)
  (vector-map (lambda (v) (vector-ref v j)) A))

(define (matrix-set! A i j val)
  (vector-set! (vector-ref A i) j val))

(define (matrix-transport A)
  (let ([n (vector-length (vector-ref A 0))])
    (vector-map (lambda (j) (matrix-ref-column A j))
                (vector-range n))))
(define (matrix-transport/parallel A)
  (let ([n (vector-length (vector-ref A 0))])
    (vector-map/parallel (lambda (j) (matrix-ref-column A j))
                         (vector-range n))))

(define (matrix-copy A)
  (let [(a (make-vector (vector-length A)))]
    (vector-for-each (lambda (j) (vector-set! a j (vector-copy (vector-ref A j))))
                     (vector-range (vector-length A)))
    a))
(define (matrix-copy/parallel A)
  (let [(a (make-vector (vector-length A)))]
    (vector-for-each/parallel (lambda (j) (vector-set! a j (vector-copy (vector-ref A j))))
                              (vector-range (vector-length A)))
    a))


;; A + B + C + ..
(define (matrix-add A . B)
  (if (null? B)
      A
      (apply vector-map
             (cons (lambda (v . vs)
                     (apply vector-map (cons + (cons v vs))))
                   (cons A B)))))
(define (matrix-add/parallel A . B)
  (if (null? B)
      A
      (apply vector-map/parallel
             (cons (lambda (v . vs)
                     (apply vector-map (cons + (cons v vs))))
                   (cons A B)))))

;; A - B - C - ..
(define (matrix-minus A . B)
  (if (null? B)
      A
      (apply vector-map
             (cons (lambda (v . vs)
                     (apply vector-map (cons - (cons v vs))))
                   (cons A B)))))

;; Ab
(define (matrix-*-vector A b)
  (vector-map (lambda (v)
                (vector-apply + (vector-map * v b)))
              A))
(define (matrix-*-vector/parallel A b)
  (vector-map/parallel (lambda (v)
                         (vector-apply + (vector-map * v b)))
                       A))

;; cA
(define (matrix-scalar-product c A)
  (vector-map (lambda (v) (vector-scalar-product c v)) A))
(define (matrix-scalar-product/parallel c A)
  (vector-map/parallel (lambda (v) (vector-scalar-product c v)) A))

;; x y^T i.e. a column vector product a row vector
(define (column-*-row x y)
  (vector-map (lambda (ele-x)
                (vector-scalar-product ele-x y))
              x))
(define (column-*-row/parallel x y)
  (vector-map/parallel (lambda (ele-x)
                (vector-scalar-product ele-x y))
                       x))

;; A B
(define (matrix-*-matrix A B)
  (vector-map (lambda (v-of-A)
                (vector-map
                 (lambda (j)
                   (vector-sum (vector-map *
                                           v-of-A
                                           (matrix-ref-column B j))))
                 (vector-range 0 (vector-length (vector-ref B 0)) 1)))
              A))
(define (matrix-*-matrix/parallel A B)
  (vector-map/parallel (lambda (v-of-A)
                         (vector-map
                          (lambda (j)
                            (vector-sum (vector-map *
                                                    v-of-A
                                                    (matrix-ref-column B j))))
                          (vector-range 0 (vector-length (vector-ref B 0)) 1)))
                       A))

